function [betas, betasRec, betasExp ] = regressions_Backus( yields, mats, rec )

% yields: T * numYields where yields are annually and continuously compounded
% mats: 1*numYields with maturities in months
% rec: 0 1 regression dummy, where 0 is boom and 1 is recression

h=mats(2)-mats(1); %horizon is equal to step between maturities

%%%%%%%%%%%%%%%%%%%%%%%%PREPARE VARIABLES%%%%%%%%%%%%%%%

T=size(yields,1);
N=size(yields,2);

matsYears = mats/12; %mats is in months
%forwards: f[t]_(n-h)_n = log(P[t]_(n-h)_n) - log(P[t]_n_n)
%                     = -(n-h)*y[t]_(n-h) - (-(n)*y[t]_(n))
%                     = n*y[t]_(n)        -(n-h)*y[t]_(n-h)
forwards  = repmat(matsYears(1,2:end),T,1).*yields(:,2:end)-repmat(matsYears(1,1:end-1),T,1).*yields(:,1:end-1);
forwards=[yields(:,1) forwards]; %f(0,h)=y(h);

%%%%%%%%%%%%%EXPECTATION HYPOTHESIS REGRESSIONS%%%%%%%%%

betas=NaN(N,2);
betasRec=NaN(N,2);
betasExp=NaN(N,2);

for i=2:N;
    %find first observation for current maturity
    sumNans = sum(isnan(yields(:,i)));
    cumSumNans = cumsum(isnan(yields(:,i)));
    firstObs = find((sumNans - [0; cumSumNans(1:end-1)])==0,1,'first');
    
    %run Backus et al. forward rate regressions
    
    y=forwards(firstObs+h:end,i-1)-yields(firstObs:end-h,1);
    x=[ones(T-h-(firstObs-1),1) (forwards(firstObs:end-h,i)-yields(firstObs:end-h,1))];
    
    [bv,sebv,R2v,R2vadj]=robustOLS(y,x,h,1); %h lags, 1 is Newey-West se
    
    betas(i,:)=[bv(2) sebv(2)];
    
    %with recession dummies
    xD=[x.*(rec(firstObs:end-h,1)*ones(1,2)) x.*(1-(rec(firstObs:end-h,1)*ones(1,2)))];
    
    [bv,sebv,R2v,R2vadj]=robustOLS(y,xD,h,1); %h lags, 1 is Newey-West se
    
    betasRec(i,:)=[bv(2) sebv(2)];
    betasExp(i,:)=[bv(4) sebv(4)];
    
end

end